Adaptive Control is Reversed Between Hands After Left Hemisphere Stroke and Lost Following Right Hemisphere Stroke¶

*Rini Varghese, James E Gordon, Robert L Sainburg, Carolee J Winstein, Nicolas Schweighofer*

Now published in the Proceedings of the National Academy of Sciences Analysis based on reviewer comments in NB section #9

In [2]:
library(IRdisplay)
display_html(
'<script>  
code_show=true; 
function code_toggle() {
  if (code_show){
    $(\'div.input\').hide();
  } else {
    $(\'div.input\').show();
  }
  code_show = !code_show
}  
$( document ).ready(code_toggle);
</script>
  <form action="javascript:code_toggle()">
    <input type="submit" value="Click here to toggle on/off the raw code.">
 </form>'
)

Load libraries & set theme¶

In [37]:
rm(list=ls())
options(warn=-1)
ReqdLibs = c("repr","ggplot2","ggpubr","patchwork","IRdisplay","dplyr","tidyr","lemon","lme4","gghalves","robustlmm",
             "lme4","ggExtra","grid","nlme","sjPlot","sjstats","table1","emmeans","broom","svglite","influence.ME","ggplotify")

invisible(lapply(ReqdLibs, library, character.only = TRUE))

thm = theme(
          strip.text.x=element_text(size=20,face="bold"),
          strip.text.y=element_text(size=20,face="bold"),
          legend.text=element_text(size=16,face="bold"),
          legend.position = "top",
          legend.title=element_text(size=16,face="bold"),
          title =element_text(size=14, face='bold'),
          text = element_text(colour = "black",size=18), 
          plot.title = element_text(colour = "black",size = 22, face = "bold"),
          axis.ticks.length = unit(0.3,"cm"),
          axis.line = element_line(colour = "black",size=0.85),
          axis.ticks = element_line(colour = "black",size=0.85),
          axis.text = element_text(colour = "black",size=24),
          axis.title=element_text(size=25)) 

display_markdown("Libraries loaded. Theme set")
options(warn=0)

Libraries loaded. Theme set

In [38]:
annotate_npc <- function(label, x, y, ...)
{
  ggplot2::annotation_custom(grid::textGrob(
    x = unit(x, "npc"), y = unit(y, "npc"), label = label, ...))
}

Hypothesis & predictions¶

Null Hypothesis: Responsibility assignment is fully flexible
Alt. Hypothesis: Responsibility assignment relies on right hemisphere specialization

In [39]:
Hand = c(rep("Left Hand", 3), rep("Right Hand", 3))
Group = rep(c("Left Hemiparesis", "Controls", "Right Hemiparesis"), 2)
Pred = c(rep("alt", 6), rep("null", 6))
#     LLHP1 LCTRL1 LRHP1 RLHP1 RCTRL1 RRHP1 LLHP0 LCTRL0 LRHP0 RLHP0 RCTRL0 RRHP0
g = c(0, -0.5, 0, 0, 0, -0.5, -0.75, -0.5, 0, 0, 0, -0.5)
acf = c(0, 0.5, 0, 0, 0, 0.5, 0.75, 0.5, 0, 0, 0, 0.5)


hypoDat = data.frame(Hand, Group, Pred, g, acf)
hypoDat$Hand = factor(hypoDat$Hand, levels = c("Left Hand", "Right Hand"))
hypoDat$Group = factor(hypoDat$Group, levels = c("Left Hemiparesis", "Controls","Right Hemiparesis"))
hypoDat
A data.frame: 12 × 5
HandGroupPredgacf
<fct><fct><chr><dbl><dbl>
Left Hand Left Hemiparesis alt 0.000.00
Left Hand Controls alt -0.500.50
Left Hand Right Hemiparesisalt 0.000.00
Right HandLeft Hemiparesis alt 0.000.00
Right HandControls alt 0.000.00
Right HandRight Hemiparesisalt -0.500.50
Left Hand Left Hemiparesis null-0.750.75
Left Hand Controls null-0.500.50
Left Hand Right Hemiparesisnull 0.000.00
Right HandLeft Hemiparesis null 0.000.00
Right HandControls null 0.000.00
Right HandRight Hemiparesisnull-0.500.50
A data.frame: 12 × 5
HandGroupPredgacf
<fct><fct><chr><dbl><dbl>
Left Hand Left Hemiparesis alt 0.000.00
Left Hand Controls alt -0.500.50
Left Hand Right Hemiparesisalt 0.000.00
Right HandLeft Hemiparesis alt 0.000.00
Right HandControls alt 0.000.00
Right HandRight Hemiparesisalt -0.500.50
Left Hand Left Hemiparesis null-0.750.75
Left Hand Controls null-0.500.50
Left Hand Right Hemiparesisnull 0.000.00
Right HandLeft Hemiparesis null 0.000.00
Right HandControls null 0.000.00
Right HandRight Hemiparesisnull-0.500.50
In [40]:
options(repr.plot.width = 12, repr.plot.height = 7)
options(warn = -1)

gr1 = ggplot(data = hypoDat, aes(x = Group, y = g-0.2, shape = Group, col = Hand, lty = Pred)) +
    stat_summary(mapping = NULL, data = NULL, show.legend = FALSE, geom = "pointrange",
        fun.data = mean_se, position = position_dodge(width = 0), size = 1.5, lwd = 3,
        alpha = 0.7) + geom_line(aes(group = interaction(Hand, Pred), col = Hand,
    lty = Pred), size = 2, alpha = 1, show.legend = FALSE, position = position_dodge(width = 0.25),
    arrow = arrow(ends = "both", type = "open")) + geom_hline(yintercept = 0, linetype = "dashed",
    size = 1) + coord_capped_cart(ylim = c(-1.38, 1.38), bottom = "both", left = "both") +
    xlab("") + ylab("gain") + scale_color_manual(values = c("#110176", "#E93324")) +
    scale_fill_manual(values = c("#110176", "#E93324")) + scale_shape_manual(values = c(15,
    16, 17)) + theme_pubclean() + thm + annotate_npc("Loss", x = unit(c(0.35), "npc"),
    y = unit(c(0.4), "npc"), gp = gpar(fontsize = 20, fontface = 4, col = "black")) +
    annotate_npc("Amplify", x = unit(c(0.14), "npc"), y = unit(c(0.25), "npc"), gp = gpar(fontsize = 20,
        fontface = 4, col = "black")) + annotate_npc("Switch", x = unit(c(0.65),
    "npc"), y = unit(c(0.25), "npc"), gp = gpar(fontsize = 20, fontface = 4, col = "black")) +

theme(axis.line.x = element_line(colour = "white"), axis.text.x = element_blank(),
    axis.ticks.x = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(),
    axis.line.y = element_line(arrow = grid::arrow(length = unit(0.7, "cm"), ends = "both")))

gr2 = ggplot(data = hypoDat, aes(x = Group, y = acf +0.2, shape = Group, col = Hand, linetype = Pred)) +
    stat_summary(mapping = NULL, data = NULL, show.legend = FALSE, geom = "pointrange",
        fun.data = mean_se, position = position_dodge(width = 0), size = 1.5, lwd = 3,
        alpha = 0.7) + geom_line(aes(group = interaction(Hand, Pred), col = Hand,
    linetype = Pred), size = 2, alpha = 1, show.legend = FALSE, position = position_dodge(width = 0.25),
    arrow = arrow(ends = "both", type = "open")) + geom_hline(yintercept = 0, linetype = "dashed",
    size = 1) + coord_capped_cart(ylim = c(-1.38, 1.38), bottom = "both", left = "both") +
    xlab("") + ylab("acf") + scale_color_manual(values = c("#110176", "#E93324")) +
    scale_fill_manual(values = c("#110176", "#E93324")) + scale_shape_manual(values = c(15,
    16, 17)) + annotate_npc("Loss", x = unit(c(0.35), "npc"), y = unit(c(0.6), "npc"),
    gp = gpar(fontsize = 20, fontface = 4, col = "black")) + annotate_npc("Amplify",
    x = unit(c(0.14), "npc"), y = unit(c(0.75), "npc"), gp = gpar(fontsize = 20,
        fontface = 4, col = "black")) + annotate_npc("Switch", x = unit(c(0.65),
    "npc"), y = unit(c(0.77), "npc"), gp = gpar(fontsize = 20, fontface = 4, col = "black")) +
    theme_pubclean() + thm + theme(axis.line.x = element_line(colour = "white"),
    axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.text.y = element_blank(),
    axis.ticks.y = element_blank(), axis.line.y = element_line(arrow = grid::arrow(length = unit(0.7,
        "cm"), ends = "both")))


gr1 + gr2 

options(warn = 0)

Participant characteristics¶

In [41]:
demo = read.csv("forTable1_v2.csv", stringsAsFactors = FALSE)
library("table1")

demo$Sex = factor(demo$Sex, levels = c("M", "F"), labels = c("Male", "Female"))
# head(demo)

rndr <- function(x, name, ...) {
    if (!is.numeric(x))
        return(render.categorical.default(x))
    what <- switch(name, Age = "Mean (SD)", EHI = "Median [Min, Max]", Chronicity = "Median [Min, Max]",
        UEFM = "Median [Min, Max]")
    parse.abbrev.render.code(c("", what))(x)
}
strata <- c(split(demo, demo$SOL))
labels <- list(variables = list(Age = "Age (years)", Sex = "Sex", EHI = "Edinburgh Handedness Inventory",
    Chronicity = "Chronicity (months)", UEFM = "UE Fugl-Meyer (/66)"))

tt2 = table1(strata, labels, render = rndr)
display_html(tt2)
Controls
(N=20)
Left Hemiparesis
(N=15)
Right Hemiparesis
(N=15)
Age (years)
Mean (SD) 26.9 (3.99) 59.4 (14.3) 63.5 (10.0)
Sex
Male 9 (45.0%) 8 (53.3%) 11 (73.3%)
Female 11 (55.0%) 7 (46.7%) 4 (26.7%)
Edinburgh Handedness Inventory
Median [Min, Max] 89.4 [44.0, 100] 100 [50.0, 100] 100 [60.0, 100]
Chronicity (months)
Median [Min, Max] NA 93.0 [23.0, 211] 95.0 [19.0, 141]
UE Fugl-Meyer (/66)
Median [Min, Max] NA 52.0 [13.0, 64.0] 55.0 [20.0, 63.0]

Within-trial analysis (g)¶

Read data¶

In [42]:
data = read.csv("YNC_AllData_120321_full.csv", stringsAsFactors = FALSE)  #

data %>%
    filter(CondG == c("BmC")) %>%
    mutate(Adj = PVDE - IDE) %>%
{
    bmCData <<- .
}

head(bmCData)
A data.frame: 6 × 33
GrpSubjSubNoCondCondGDistTrialHandCTRT⋯FDISTCURVOnAxErrOffAxErrXPosPVYPosPVXPosOffYPosOffXDistPVAdj
<chr><chr><int><chr><chr><int><int><chr><dbl><dbl>⋯<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
1YNCSubj0011BmC05BmC51left0.919-0.03107⋯6.280590.15171 3.162355 0.11256510.92261551.0414480.92311461.073100-0.004901704-5.20114
2YNCSubj0011BmC05BmC52left0.982-0.10269⋯4.941050.08149 1.562621 1.07743700.93467581.0260620.93276331.057103 0.007158609-3.49286
3YNCSubj0011BmC05BmC53left0.928-0.03145⋯7.036780.16236 1.992743 1.17369290.92840341.0326990.93372591.061404 0.000886155-5.48651
4YNCSubj0011BmC05BmC54left0.992-0.10407⋯6.597750.28700 2.478371-1.51805080.91091581.0357500.90680851.066261-0.016601443-4.83403
5YNCSubj0011BmC05BmC55left0.946-0.05807⋯3.976410.29868-0.291631-0.20045630.91089481.0210320.91998441.038560-0.016622415-0.18322
6YNCSubj0011BmC05BmC56left0.994-0.09745⋯5.602930.32918 1.634333-2.56340310.90296151.0288110.89635491.057820-0.024555732-4.46440

Outlier detection & removal¶

In [43]:
isnt_out_z <- function(x, thres = 3, na.rm = TRUE) {
  abs(x - mean(x, na.rm = na.rm)) <= thres * sd(x, na.rm = na.rm)
}

compute_group_non_outliers <- . %>%
  # Compute per group mean values of columns
  group_by(group) %>%
  summarise_if(is.numeric, mean) %>%
  ungroup() %>%

  # Detect outliers among groups
  mutate_if(is.numeric, isnt_out_z)  %>% 
  # Remove unnecessary columns
  select_if(Negate(is.numeric))
In [44]:
head(bmCData)
colnames(bmCData)[1]="Grp"
bmCData  %>% 
unite(col = "group", Grp,Subj,Cond,Trial) %>% 
compute_group_non_outliers() %>% 
filter(Adj == FALSE | PV==FALSE | IDE==FALSE & PRDE==FALSE | PVDE==FALSE)  %>%  #
{.$group->> outs}
cat(paste(shQuote(outs, type="cmd"), collapse=", "))

dfTrials = dim(bmCData)[1]/3
perc_rem = (length(outs)/dfTrials)*100 # we divide by 3 because the data frame consists of 3 values (for hands & cursor) per trial

display_markdown(paste("<br>**With these outliers, we will be removing",length(outs),"data points out of ",dfTrials,"(",round(perc_rem,3),"%) of the data**"))


#below printout to rem
A data.frame: 6 × 33
GrpSubjSubNoCondCondGDistTrialHandCTRT⋯FDISTCURVOnAxErrOffAxErrXPosPVYPosPVXPosOffYPosOffXDistPVAdj
<chr><chr><int><chr><chr><int><int><chr><dbl><dbl>⋯<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
1YNCSubj0011BmC05BmC51left0.919-0.03107⋯6.280590.15171 3.162355 0.11256510.92261551.0414480.92311461.073100-0.004901704-5.20114
2YNCSubj0011BmC05BmC52left0.982-0.10269⋯4.941050.08149 1.562621 1.07743700.93467581.0260620.93276331.057103 0.007158609-3.49286
3YNCSubj0011BmC05BmC53left0.928-0.03145⋯7.036780.16236 1.992743 1.17369290.92840341.0326990.93372591.061404 0.000886155-5.48651
4YNCSubj0011BmC05BmC54left0.992-0.10407⋯6.597750.28700 2.478371-1.51805080.91091581.0357500.90680851.066261-0.016601443-4.83403
5YNCSubj0011BmC05BmC55left0.946-0.05807⋯3.976410.29868-0.291631-0.20045630.91089481.0210320.91998441.038560-0.016622415-0.18322
6YNCSubj0011BmC05BmC56left0.994-0.09745⋯5.602930.32918 1.634333-2.56340310.90296151.0288110.89635491.057820-0.024555732-4.46440
"LHP_Subj003_BmC05_17", "LHP_Subj004_BmC15_6", "LHP_Subj005_BmC15_17", "LHP_Subj007_BmC05_1", "LHP_Subj007_BmC05_4", "LHP_Subj009_BmC05_10", "LHP_Subj009_BmC05_18", "LHP_Subj009_BmC05_5", "LHP_Subj009_BmC05_8", "LHP_Subj009_BmC15_11", "LHP_Subj011_BmC05_10", "LHP_Subj011_BmC05_12", "LHP_Subj011_BmC05_13", "LHP_Subj011_BmC05_15", "LHP_Subj011_BmC05_17", "LHP_Subj011_BmC05_18", "LHP_Subj011_BmC05_19", "LHP_Subj011_BmC05_20", "LHP_Subj011_BmC15_1", "LHP_Subj011_BmC15_11", "LHP_Subj011_BmC15_12", "LHP_Subj011_BmC15_15", "LHP_Subj011_BmC15_16", "LHP_Subj011_BmC15_19", "LHP_Subj011_BmC15_2", "LHP_Subj012_BmC05_9", "LHP_Subj012_BmC15_7", "LHP_Subj013_BmC05_1", "LHP_Subj013_BmC05_10", "LHP_Subj013_BmC05_11", "LHP_Subj013_BmC05_14", "LHP_Subj013_BmC05_15", "LHP_Subj013_BmC05_19", "LHP_Subj013_BmC05_3", "LHP_Subj013_BmC05_5", "LHP_Subj013_BmC05_7", "LHP_Subj013_BmC05_8", "LHP_Subj013_BmC05_9", "LHP_Subj013_BmC15_20", "LHP_Subj014_BmC05_15", "LHP_Subj016_BmC15_8", "RHP_Subj001_BmC05_13", "RHP_Subj001_BmC15_2", "RHP_Subj002_BmC05_7", "RHP_Subj005_BmC15_15", "RHP_Subj006_BmC15_14", "RHP_Subj007_BmC05_2", "RHP_Subj007_BmC15_6", "RHP_Subj008_BmC05_15", "RHP_Subj011_BmC05_3", "RHP_Subj011_BmC05_8", "RHP_Subj013_BmC05_12", "RHP_Subj013_BmC05_13", "RHP_Subj013_BmC05_8", "YNC_Subj004_BmC05_16"


With these outliers, we will be removing 55 data points out of 2000 ( 2.75 %) of the data

In [45]:
# rem2 = c("LHP_Subj002_BmC05_10", "LHP_Subj003_BmC05_17", "LHP_Subj004_BmC05_5", "LHP_Subj004_BmC15_20", "LHP_Subj004_BmC15_6", "LHP_Subj005_BmC15_17", "LHP_Subj007_BmC05_1", "LHP_Subj007_BmC05_4", "LHP_Subj007_BmC15_18", "LHP_Subj008_BmC05_10", "LHP_Subj009_BmC05_10", "LHP_Subj009_BmC05_13", "LHP_Subj009_BmC05_18", "LHP_Subj009_BmC05_4", "LHP_Subj009_BmC05_5", "LHP_Subj009_BmC05_8", "LHP_Subj009_BmC15_10", "LHP_Subj009_BmC15_11", "LHP_Subj009_BmC15_16", "LHP_Subj011_BmC05_10", "LHP_Subj011_BmC05_17", "LHP_Subj011_BmC05_18", "LHP_Subj011_BmC05_19", "LHP_Subj011_BmC05_20", "LHP_Subj011_BmC05_6", "LHP_Subj011_BmC05_9", "LHP_Subj011_BmC15_1", "LHP_Subj011_BmC15_11", "LHP_Subj011_BmC15_12", "LHP_Subj011_BmC15_13", "LHP_Subj011_BmC15_15", "LHP_Subj011_BmC15_16", "LHP_Subj011_BmC15_19", "LHP_Subj011_BmC15_2", "LHP_Subj012_BmC05_6", "LHP_Subj012_BmC05_9", "LHP_Subj012_BmC15_7", "LHP_Subj012_BmC15_8",  "LHP_Subj013_BmC05_10",  "LHP_Subj013_BmC05_9", "LHP_Subj013_BmC05_15", "LHP_Subj013_BmC05_19", "LHP_Subj013_BmC05_3", "LHP_Subj013_BmC05_5", "LHP_Subj013_BmC05_7", "LHP_Subj013_BmC05_8",  "LHP_Subj014_BmC05_15", "LHP_Subj014_BmC05_6", "LHP_Subj016_BmC15_8", "RHP_Subj001_BmC05_13", "RHP_Subj001_BmC15_14", "RHP_Subj001_BmC15_2", "RHP_Subj002_BmC05_7", "RHP_Subj003_BmC05_17", "RHP_Subj003_BmC05_8", "RHP_Subj003_BmC15_8", "RHP_Subj005_BmC15_15", "RHP_Subj006_BmC05_3", "RHP_Subj006_BmC05_8", "RHP_Subj006_BmC15_14", "RHP_Subj006_BmC15_7", "RHP_Subj007_BmC05_2", "RHP_Subj007_BmC15_3", "RHP_Subj007_BmC15_6", "RHP_Subj008_BmC05_15", "RHP_Subj011_BmC05_16", "RHP_Subj011_BmC05_3", "RHP_Subj011_BmC05_8", "RHP_Subj012_BmC15_15", "RHP_Subj012_BmC15_6", "RHP_Subj013_BmC05_12", "RHP_Subj013_BmC05_13", "RHP_Subj013_BmC05_17", "RHP_Subj013_BmC05_8", "RHP_Subj015_BmC05_5", "YNC_Subj004_BmC05_16", "YNC_Subj009_BmC05_4"
# )

rem2 = c("LHP_Subj003_BmC05_17", "LHP_Subj004_BmC15_6", "LHP_Subj005_BmC15_17", "LHP_Subj007_BmC05_1", "LHP_Subj007_BmC05_4", "LHP_Subj009_BmC05_10", "LHP_Subj009_BmC05_18", "LHP_Subj009_BmC05_5", "LHP_Subj009_BmC05_8", "LHP_Subj009_BmC15_11", "LHP_Subj011_BmC05_10", "LHP_Subj011_BmC05_12", "LHP_Subj011_BmC05_13", "LHP_Subj011_BmC05_15", "LHP_Subj011_BmC05_17", "LHP_Subj011_BmC05_18", "LHP_Subj011_BmC05_19", "LHP_Subj011_BmC05_20", "LHP_Subj011_BmC15_1", "LHP_Subj011_BmC15_11", "LHP_Subj011_BmC15_12", "LHP_Subj011_BmC15_15", "LHP_Subj011_BmC15_16", "LHP_Subj011_BmC15_19", "LHP_Subj011_BmC15_2", "LHP_Subj012_BmC05_9", "LHP_Subj012_BmC15_7", "LHP_Subj013_BmC05_1", "LHP_Subj013_BmC05_10", "LHP_Subj013_BmC05_11", "LHP_Subj013_BmC05_14", "LHP_Subj013_BmC05_15", "LHP_Subj013_BmC05_19", "LHP_Subj013_BmC05_3", "LHP_Subj013_BmC05_5", "LHP_Subj013_BmC05_7", "LHP_Subj013_BmC05_8", "LHP_Subj013_BmC05_9", "LHP_Subj013_BmC15_20", "LHP_Subj014_BmC05_15", "LHP_Subj016_BmC15_8", "RHP_Subj001_BmC05_13", "RHP_Subj001_BmC15_2", "RHP_Subj002_BmC05_7", "RHP_Subj005_BmC15_15", "RHP_Subj006_BmC15_14", "RHP_Subj007_BmC05_2", "RHP_Subj007_BmC15_6", "RHP_Subj008_BmC05_15", "RHP_Subj011_BmC05_3", "RHP_Subj011_BmC05_8", "RHP_Subj013_BmC05_12", "RHP_Subj013_BmC05_13", "RHP_Subj013_BmC05_8", "YNC_Subj004_BmC05_16"
)
In [46]:
bmCData  %>% 
unite(col = "group",Grp,Subj,Cond,Trial) %>% 
filter(!group %in% rem2) %>% 
separate(col = group, into = c("Grp","Subj","Cond","Trial"),sep = "_")  %>% 
pivot_wider(id_cols = NULL, names_from= Hand,values_from = c(CT,RT,PPT,TPV,MT,PV,EE,VE,IDE,PRDE,PVDE,FDE,FPE,IDIST,PVDIST,FDIST,
                                             CURV,OnAxErr,OffAxErr,Adj,XPosPV,YPosPV,XPosOff,YPosOff,XDistPV))  %>% 

pivot_longer(cols = c(Adj_left, Adj_right), names_to = "Hand",values_to = "Adjustment") %>% 
select(Grp,Subj,Cond,Trial,SubNo,CondG,Dist,Hand,IDE_cursor,Adjustment)  %>% 

{.->>bmCData2}

bmCData2$Hand <- factor(bmCData2$Hand, levels = c("Adj_left", "Adj_right"), labels = c("Left", "Right")) 

head(bmCData2)
bmCData2$Grp = factor(bmCData2$Grp, levels=c('YNC','RHP','LHP'))
A tibble: 6 × 10
GrpSubjCondTrialSubNoCondGDistHandIDE_cursorAdjustment
<chr><chr><chr><chr><int><chr><int><fct><dbl><dbl>
YNCSubj001BmC0511BmC5Left -7.20574-5.20114
YNCSubj001BmC0511BmC5Right-7.20574 6.50602
YNCSubj001BmC0521BmC5Left -6.89997-3.49286
YNCSubj001BmC0521BmC5Right-6.89997 2.20555
YNCSubj001BmC0531BmC5Left -2.80350-5.48651
YNCSubj001BmC0531BmC5Right-2.80350 6.49121

g Single Represesentative Subjects¶

In [47]:
options(repr.plot.width = 13.5, repr.plot.height = 4.5)

Ind_subjGL = 
ggplot(bmCData2 %>% filter(SubNo==36), aes(y = IDE_cursor, x = Adjustment,na.rm = TRUE,group = Hand)) +
       geom_point(aes(col = Hand),stroke = 1.5, size = 3.5, shape = 15, show.legend = FALSE,alpha = 0.3) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'****'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),label.x = -20, label.y = c(20,14), size=8,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-40,20),bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand),  method = 'lm', fill=NA,formula = y~x, lwd=3.5,show.legend = FALSE,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + 
       scale_alpha_manual(values = c(0.7, 0.2)) +
       theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjL.svg', plot=Ind_subjGL, width=4.5, height=4.5)

Ind_subjG = 
ggplot(bmCData2 %>% filter(SubNo==14), aes(y = Adjustment, x = IDE_cursor,na.rm = TRUE)) +
       geom_point(aes(col = Hand),stroke = 1.5, alpha = 0.3, size = 3.5, shape = 16, show.legend = FALSE) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'***'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),label.x = -20, label.y = c(20,17),size=7.5,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-12,20), bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand,fill=Hand), fill=NA,method = 'lm', formula = y~x, lwd=3.5,show.legend = FALSE,alpha=0.1,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjG.svg', plot=Ind_subjG, width=4.5, height=4.5)

Ind_subjGR = 
ggplot(bmCData2 %>% filter(SubNo==21), aes(y = IDE_cursor, x =Adjustment ,na.rm = TRUE,group = Hand)) +
       geom_point(aes(col = Hand),stroke = 1.5, size = 3.5, shape = 17, show.legend = FALSE,alpha = 0.3) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'***'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),label.x = -1, label.y = c(20,17), size=7.5,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-10,20),bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand),  method = 'lm', fill=NA,formula = y~x, lwd=3.5,show.legend = FALSE,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + 
       theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjR.svg', plot=Ind_subjGR, width=4.5, height=4.5)

Ind_subjGALL = Ind_subjGL | Ind_subjG | Ind_subjGR

# ggsave(file='Ind_subjGALL.svg', plot=Ind_subjGALL, width=13.5, height=4.5)
Ind_subjGALL
In [48]:
options(dplyr.summarise.inform = FALSE)
options(warn = -1)
bmCData2 %>%
    group_by(Grp) %>%
    summarise(mean_theta_cursor = mean(IDE_cursor, na.rm = TRUE), se_acf = se(IDE_cursor, na.rm = TRUE))
options(dplyr.summarise.inform = TRUE)
A tibble: 3 × 3
Grpmean_theta_cursorse_acf
<fct><dbl><dbl>
YNC-4.4953510.1975443
RHP-1.8891870.3112893
LHP-4.5862250.4026761

Group-wise estimates for g (from fixed effects)¶

In [49]:
options(warn = -1)
mixedModEstInter = suppressMessages(lmer(data = bmCData2, IDE_cursor ~ Grp + Grp:Hand:Adjustment +
    (1 + Grp:Hand:Adjustment | SubNo)))
coefs <- coef(summary(mixedModEstInter))
data.frame(coefs[, 1:2])
options(warn = 0)
A data.frame: 9 × 2
EstimateStd..Error
<dbl><dbl>
(Intercept)-4.62643231.30975520
GrpRHP 2.93583071.99869422
GrpLHP 0.40704872.00068142
GrpYNC:HandLeft:Adjustment-0.54746030.06922628
GrpRHP:HandLeft:Adjustment-0.11555440.08116438
GrpLHP:HandLeft:Adjustment-0.15357230.06743758
GrpYNC:HandRight:Adjustment-0.14423720.07575751
GrpRHP:HandRight:Adjustment-0.24884560.04458200
GrpLHP:HandRight:Adjustment-0.18757440.10694754

Subject-wise estimates for g (from random effects)¶

In [50]:
Left = c(coef(mixedModEstInter)$SubNo[1:20,"GrpYNC:HandLeft:Adjustment"],
         coef(mixedModEstInter)$SubNo[21:32,"GrpRHP:HandLeft:Adjustment"],
         coef(mixedModEstInter)$SubNo[44:47,"GrpRHP:HandLeft:Adjustment"],         
         coef(mixedModEstInter)$SubNo[33:43,"GrpLHP:HandLeft:Adjustment"],
         coef(mixedModEstInter)$SubNo[48:50,"GrpLHP:HandLeft:Adjustment"])

Right = c(coef(mixedModEstInter)$SubNo[1:20,"GrpYNC:HandRight:Adjustment"],
         coef(mixedModEstInter)$SubNo[21:32,"GrpRHP:HandRight:Adjustment"],
         coef(mixedModEstInter)$SubNo[44:47,"GrpRHP:HandRight:Adjustment"],         
         coef(mixedModEstInter)$SubNo[33:43,"GrpLHP:HandRight:Adjustment"],
         coef(mixedModEstInter)$SubNo[48:50,"GrpLHP:HandRight:Adjustment"])


Grp = c(matrix('YNC',20),matrix('RHP',15),matrix('LHP',15)) #,

Subj = c(1:length(Left))
Slopes = data.frame(Subj,Grp,Left,Right)

# Pivot data table for analysis & plotting
Slopes  %>% 
pivot_longer(cols = c(Left,Right),names_to = "Hand",values_to = "g")  %>% 
{.->>MMEstGains}
MMEstGains$Grp = factor(MMEstGains$Grp, levels=c('LHP','YNC','RHP'))
head(MMEstGains)
A tibble: 6 × 4
SubjGrpHandg
<int><fct><chr><dbl>
1YNCLeft -0.54379971
1YNCRight-0.06728542
2YNCLeft -0.62510505
2YNCRight-0.12980678
3YNCLeft -0.37560941
3YNCRight 0.14613980

Plotting g (hypothesis & results)¶

In [51]:
options(repr.plot.width = 11, repr.plot.height = 7)
options(warn = -1)

rs1 = ggplot(data = MMEstGains, aes(x = Grp, y = g, shape = Grp, col = Hand)) + stat_summary(mapping = NULL,
    data = NULL, fun.data = mean_se, alpha = 0.8, geom = "pointrange", position = position_dodge(width = 0.5),
    size = 1.5, lwd = 3, show.legend = FALSE) + geom_half_point_panel(aes(col = Hand),
    size = 3, alpha = 0.3, position = position_dodge(width = 0.5), show.legend = FALSE) +
    geom_line(aes(group = Hand, col = Hand), size = 2, fun.data = mean_se, alpha = 1,
        position = position_dodge(width = 0.5), stat = "summary", show.legend = FALSE) +
    geom_hline(yintercept = 0, linetype = "dashed", size = 1) + coord_capped_cart(ylim = c(-0.75,
    0.75), bottom = "both", left = "both") + xlab("") + ylab("g") + annotate_npc("paretic left",
    x = unit(c(0.08), "npc"), y = unit(c(0.3), "npc"), rot = 90, gp = gpar(fontsize = 20,
        fontface = 3, col = "#110176")) + annotate_npc("paretic right", x = unit(c(0.85),
    "npc"), y = unit(c(0.2), "npc"), rot = 90, gp = gpar(fontsize = 20, fontface = 3,
    col = "#E93324")) + scale_color_manual(values = c("#110176", "#E93324")) + scale_fill_manual(values = c("#110176",
    "#E93324")) + scale_shape_manual(values = c(15, 16, 17)) + theme_pubclean() +
    thm + theme(axis.line.x = element_line(colour = "white"), axis.text.x = element_blank(),
    axis.ticks.x = element_blank(), legend.position = "none")


Fig1 = gr1 + rs1

# ggsave(file = "Fig1_g.svg", plot = Fig1, width = 11, height = 7)

Fig1
options(warn = 0)

Simple LMER for g & Model Diagnostics¶

In [52]:
gMod = lmer(g ~ 1 + Hand + Grp + Hand*Grp + (1| Subj), MMEstGains)
estex = influence(gMod,"Subj")

options(repr.plot.width = 14, repr.plot.height = 7)

diagGp1 = as.ggplot(~plot(estex, which="cook",cutoff=.04, sort=TRUE,xlab="Cook's Distance",ylab = "Subj"))
diagGp2 = as.ggplot(~plot(estex, which="dfbetas",cutoff=.28, xlab="DFBetas",ylab = "Subj"))

diagGp1|diagGp2
In [53]:
estex.obs <- influence(gMod, obs=TRUE) 
cks.d <- cooks.distance(estex.obs, parameter=3) 
outliers = which(cks.d > 4/100)

gMod2 = exclude.influence(gMod,"Subj",outliers)

Robust LMER to compare g between groups and hands¶

Comparing g between YNC & LHP¶

In [54]:
# SET REFERENCE
MMEstGains$Hand <- factor(MMEstGains$Hand, ordered = FALSE)
MMEstGains$Hand <- relevel(MMEstGains$Hand, "Right")
MMEstGains$Grp <- relevel(MMEstGains$Grp, "LHP")

# RUN ROBUST MODEL
robust.model1L = rlmer(g ~ Hand + Grp + Hand:Grp + (1 | Subj), MMEstGains %>% filter(Grp != "RHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model1L, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model1L, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftLHP-0.025834990.03128241Inf-0.82586334.088816e-01
2Right - LeftYNC 0.402188680.02709136Inf14.84564427.423821e-50
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1LHP - YNCRight-0.032296140.05213391Inf-0.61948445.355973e-01
2LHP - YNCLeft 0.395727530.05213391Inf 7.59059833.184317e-14

Comparing g between YNC & RHP¶

In [55]:
# SET REFERENCE
MMEstGains$Hand <- factor(MMEstGains$Hand, ordered = FALSE)
MMEstGains$Hand <- relevel(MMEstGains$Hand, "Right")
MMEstGains$Grp <- relevel(MMEstGains$Grp, "RHP")

# RUN ROBUST MODEL
robust.model1R = rlmer(g ~ Hand + Grp + Hand:Grp + (1 | Subj), MMEstGains %>% filter(Grp != "LHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model1R, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model1R, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftRHP-0.16490340.04858103Inf-3.3943986.877959e-04
2Right - LeftYNC 0.38758220.04207241Inf 9.2122663.193107e-20
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1RHP - YNCRight-0.091193090.0454434Inf-2.006744.477734e-02
2RHP - YNCLeft 0.461292520.0454434Inf10.150933.282355e-24

Trial-to-trial analysis (acf)¶

Read Data¶

In [56]:
colnames(data)[1] = "Grp"
data  %>% 
filter(CondG==c("BmC")) %>% 
filter(Hand!="cursor") %>% 
mutate(Trial = as.numeric(Trial)) %>% 
select(Grp,Subj,SubNo,Cond,Dist,Hand,Trial,PVDE,XPosPV,YPosPV) %>% 

unite(col = "group",Grp,Subj,SubNo,Cond,Dist,Hand) %>% 
select(group,Trial,PVDE,XPosPV,YPosPV)   %>% 
{.->> ctrlDataF}

head(ctrlDataF)
A data.frame: 6 × 5
groupTrialPVDEXPosPVYPosPV
<chr><dbl><dbl><dbl><dbl>
1YNC_Subj001_1_BmC05_5_left1-10.324520.92261551.041448
2YNC_Subj001_1_BmC05_5_left2 -6.417640.93467581.026062
3YNC_Subj001_1_BmC05_5_left3-10.417850.92840341.032699
4YNC_Subj001_1_BmC05_5_left4-19.343860.91091581.035750
5YNC_Subj001_1_BmC05_5_left5 -9.369320.91089481.021032
6YNC_Subj001_1_BmC05_5_left6-22.259100.90296151.028811

Computing Mahalanobis distance & lagged PVDE to determine fast learning¶

In [57]:
names = unique(ctrlDataF$group)

mahals = vector()
datLags = vector()

for (i in names)

    {dat = ctrlDataF[ctrlDataF$group == i,]
     
     datLag = lag(dat$PVDE,1)
     datLags = c(datLags,datLag)
     
     mahal = mahalanobis(dat[,4:5], 
                    colMeans(dat[,4:5]),
                    cov(dat[,4:5]))
     mahals = c(mahals,mahal)
    }

lagPVDE = datLags

ctrlDataNew = data.frame(ctrlDataF,lagPVDE,mahals)

ctrlDataNew  %>% 
separate(col = group, into = c("Grp","Subj","SubNo","Cond","Dist","Hand"),sep = "_")   %>% 
{.->>ctrlDataNew4reg}

ctrlDataNew4reg$Grp = factor(ctrlDataNew4reg$Grp, levels=c('YNC','RHP','LHP'))
head(ctrlDataNew4reg)
A data.frame: 6 × 12
GrpSubjSubNoCondDistHandTrialPVDEXPosPVYPosPVlagPVDEmahals
<fct><chr><chr><chr><chr><chr><dbl><dbl><dbl><dbl><dbl><dbl>
1YNCSubj0011BmC055left1-10.324520.92261551.041448 NA6.186850
2YNCSubj0011BmC055left2 -6.417640.93467581.026062-10.324524.462720
3YNCSubj0011BmC055left3-10.417850.92840341.032699 -6.417643.160088
4YNCSubj0011BmC055left4-19.343860.91091581.035750-10.417852.501502
5YNCSubj0011BmC055left5 -9.369320.91089481.021032-19.343861.138511
6YNCSubj0011BmC055left6-22.259100.90296151.028811 -9.369322.812689

Select PVDE for ACF analysis¶

In [58]:
ctrlDataNew4reg  %>% 
filter(Trial!=1:5) %>%
select(SubNo,Trial,Dist,Hand,PVDE) %>% 
pivot_wider(id_cols=NULL,names_from = Hand,values_from = PVDE) %>% 
{.->>ctrlDataNew4ACF}

head(ctrlDataNew4ACF)
A tibble: 6 × 5
SubNoTrialDistleftright
<chr><dbl><chr><dbl><dbl>
1 65-22.25910 2.23242
1 75-17.25457 8.54909
1 85-18.71507 4.56551
1 95-23.68392 9.18148
1105-19.8991316.52577
1115 -7.35605 8.36514

All subjects ACFs¶

In [61]:
acfL = vector()
acfR = vector()

SubjID = vector()

for (i in ctrlDataNew4ACF$SubNo) {
    acfL[i] = acf(ctrlDataNew4ACF$left[ctrlDataNew4ACF$SubNo == i], lag.max = 1,
        plot = FALSE)$acf[2]
    acfR[i] = acf(ctrlDataNew4ACF$right[ctrlDataNew4ACF$SubNo == i], lag.max = 1,
        plot = FALSE)$acf[2]
    SubjID[i] = i
}


acf = data.frame(acfL, acfR)

colnames(acf)[1] <- "Left"
colnames(acf)[2] <- "Right"
# colnames(acf)[3] <- 'Cursor'


j = c(matrix("YNC", 20), matrix("RHP", 15), matrix("LHP", 15))  #,

acfData = data.frame(SubjID, data.frame(j), acf)

colnames(acfData)[2] = "Grp"

acfData %>%
    pivot_longer(cols = c("Left", "Right"), "Hand", values_to = "acf_coeff") %>%

{
    acfData2 <<- .
}

acfData2$Grp = factor(acfData2$Grp, levels = c("LHP", "YNC", "RHP"))

acf Single Represesentative Subject¶

In [62]:
eg = 1;
remT = c(1:5)

ctrlDataNew4reg %>% 
select(SubNo,Trial,Dist,Hand,PVDE,lagPVDE) %>% 
filter(SubNo==eg) %>% 

{.->> sub1Dat}
head(sub1Dat,3)


lefteg = acf(ctrlDataNew4ACF$left[ctrlDataNew4ACF$SubNo == eg], plot = FALSE)
righteg = acf(ctrlDataNew4ACF$right[ctrlDataNew4ACF$SubNo == eg], plot = FALSE)
subjACF = data.frame(lefteg$lag, lefteg$acf, righteg$acf)


colnames(subjACF)[1] <- "Lag"
colnames(subjACF)[2] <- "Left"
colnames(subjACF)[3] <- "Right"

subjACF %>%
    pivot_longer(cols = c(Left, Right), names_to = "Hand", values_to = "acf") %>% 
    filter(Lag<8) %>% 
    {subjACF <<- .}
A data.frame: 3 × 6
SubNoTrialDistHandPVDElagPVDE
<chr><dbl><chr><chr><dbl><dbl>
1115left-10.32452 NA
2125left -6.41764-10.32452
3135left-10.41785 -6.41764
In [63]:
options(repr.plot.width = 10, repr.plot.height = 5)


#--------------------ACF PLOT-----------------------#


IndFig = ggplot(data = subjACF, aes(x = Lag, y = acf, fill = Hand)) + 
geom_bar(width = 0.9, position = position_dodge(width = 0.8), stat = "identity", show.legend = FALSE,
    alpha = 0.3) + geom_bar(data = subset(subjACF, Lag == 1), aes(x = Lag, y = acf,
    fill = Hand), width = 1, position = position_dodge(width = 0.9), stat = "identity",
    show.legend = FALSE, alpha = 1) + annotate("segment", x = 0.8, xend = 0.8, y = 0.5, yend = 0.95,
           colour = "#110176", size = 2.2, arrow = arrow()) + #type = "closed",angle = 20,ends = "last",length = unit(0.06, "npc")
annotate("segment", x = 1.2, xend = 1.2, y = -0.09, yend = -0.25,
           colour = "#E93324", size = 2.2, arrow = arrow()) +
annotate_npc("lag = 1", x = unit(c(0.28), "npc"), y = unit(c(0.5), "npc"), rot = 30, gp = gpar(fontsize = 20, fontface = 3,
        col = "black")) +
coord_capped_cart(ylim = c(-0.35, 1.2), xlim = c(0, 8), bottom = "both", left = "both") + 
    xlab("lag") + ylab("acf") + theme_pubclean() + thm + scale_color_manual(values = c("#110176",
    "#E93324")) + scale_fill_manual(values = c("#110176", "#E93324"))
# ggsave(file='IndFig_acf.svg', plot=IndFig, width=5, height=5)


#-------------------- t by (t-1) plot-----------------------#


IndFig2 = ggplot(sub1Dat %>%filter(!Trial%in%remT), aes(y = PVDE, x = lagPVDE,na.rm = TRUE)) +
       geom_point(aes(col = Hand,alpha = Trial),stroke = 0.5, size = 3.5, shape = 16, show.legend = FALSE) + 
       geom_hline(yintercept = 0,size = 1,col="gray30") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'****'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),label.x = 0, label.y = c(-21,-17), size=7.4,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-25,25),bottom = 'both', left = 'both') + 
       geom_smooth(aes(col=Hand,fill=Hand),  method = 'lm', formula = y~x, lwd=2.5,show.legend = FALSE,alpha=0.1) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_alpha("Trial",range = c(0,1),breaks = c(1:20))+
       scale_fill_manual(values=c("#110176","#E93324")) + theme_pubclean() + thm +
       xlab(expression(bold(theta) [italic(peak~vel)~(t-1)])) + 
       ylab(expression(bold(theta) [italic(peak~vel)~(t)]))

# ggsave(file="IndFig2_tbyt-1.svg", plot=IndFig2, width=5, height=5)



IndFig | IndFig2

Mean ± SE for ACFs per group¶

In [64]:
options(dplyr.summarise.inform = FALSE)
options(warn = -1)
acfData2 %>%
    group_by(Grp,Hand) %>%
    summarise(mean_acf = mean(acf_coeff, na.rm = TRUE), se_acf = se(acf_coeff, na.rm = TRUE))

options(dplyr.summarise.inform = TRUE)
options(warn = 0)
A grouped_df: 6 × 4
GrpHandmean_acfse_acf
<fct><chr><dbl><dbl>
LHPLeft 0.118308900.04667883
LHPRight 0.082824710.06191646
YNCLeft 0.252820720.04616468
YNCRight 0.044513170.04740744
RHPLeft -0.014023820.05812553
RHPRight 0.210974370.04067150

Plotting acf (hypothesis & results)¶

In [65]:
options(repr.plot.width=11, repr.plot.height=7)

options(warn=-1)

rs2 = ggplot(data = acfData2,aes(x = Grp, y = acf_coeff,shape=Grp,col=Hand)) +
      stat_summary(mapping = NULL, data = NULL,show.legend = FALSE,fun.data = mean_se, alpha = 0.8,geom = "pointrange", position = position_dodge(width=0.5),size=1.5,lwd=3) +
      geom_half_point_panel(aes(col=Hand),position = position_dodge(width=0.5),size=3, alpha = 0.2,show.legend = FALSE) +
      geom_line(aes(group=Hand,col=Hand),size=2,alpha = 1,stat = "summary",
                fun.data = mean_se, position = position_dodge(width=0.5),show.legend = FALSE) + 
      geom_hline(yintercept = 0,linetype = "dashed",size = 1) + 
      coord_capped_cart(ylim = c(-0.6,0.6), bottom = 'both', left = 'both') +
      xlab("") + ylab("acf")+ 
      scale_color_manual(values=c("#110176","#E93324")) + 
      scale_fill_manual(values=c("#110176","#E93324")) + 
      scale_shape_manual(values = c(15,16,17)) + 
      annotate_npc('paretic left', x = unit(c(0.08), "npc"), y = unit(c(0.67), "npc"),rot = 90, gp = gpar(fontsize = 20, fontface = 3,col = "#110176")) +
      annotate_npc('paretic right', x = unit(c(0.93), "npc"), y = unit(c(0.73), "npc"),rot = 90, gp = gpar(fontsize = 20, fontface = 3,col = "#E93324")) +
      theme_pubclean() + thm +
      theme(axis.line.x = element_line(colour = "white"),
            axis.text.x=element_blank(),
            axis.ticks.x=element_blank(),
            legend.position = "None")

Fig2m = gr2 + rs2

# ggsave(file="Fig2_acf.svg", plot=Fig2m, width=11, height=7)

Fig2m

Simple LMER for acf & Model Diagnostics¶

In [66]:
#LMER Model
acfMod = lmer(acf_coeff ~ 1 + Hand + Grp + Hand:Grp + (1| SubjID), acfData2)

# INFLUENTIAL DATA
estex = influence(acfMod,"SubjID")
estex.obs <- influence(acfMod, obs=TRUE) 
cks.d <- cooks.distance(estex.obs, parameter=3) 
outliers = which(cks.d > (4/70))


#PLOT
options(repr.plot.width = 14, repr.plot.height = 7)
diagACFp1 = as.ggplot(~plot(estex, which="cook",cutoff=.06, sort=TRUE,xlab="Cook's Distance",ylab = "Subj"))
diagACFp2 = as.ggplot(~plot(estex, which="dfbetas",cutoff=.28, xlab="DFBetas",ylab = "Subj"))

diagACFp1|diagACFp2
In [67]:
options(repr.plot.width = 15, repr.plot.height = 5.5)

dg1 = as.ggplot(plot(acfMod))
dg2 = as.ggplot(~qqnorm(residuals(acfMod)))
dg3 = as.ggplot(~hist(residuals(acfMod)))

dg1 | dg2 | dg3
shapiro.test(residuals(acfMod))
	Shapiro-Wilk normality test

data:  residuals(acfMod)
W = 0.97471, p-value = 0.05119
In [68]:
acfMod2 = exclude.influence(acfMod,"SubjID",outliers)
# summary(acfMod2)

Robust LMER to compare acf between groups and hands¶

Comparing acf between YNC & LHP¶

In [70]:
# SET REFERENCE
acfData2$Hand <- factor(acfData2$Hand, ordered = FALSE)
acfData2$Hand <- relevel(acfData2$Hand, "Right")
acfData2$Grp <- relevel(acfData2$Grp, "LHP")

# RUN ROBUST MODEL
robust.model2L = rlmer(acf_coeff ~ Hand + Grp + Hand:Grp + (1 | SubjID), acfData2 %>%
    filter(Grp != "RHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model2L, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model2L, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftLHP-0.043907670.07320295Inf-0.59980730.5486346439
2Right - LeftYNC-0.225788070.06339561Inf-3.56157250.0003686403
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1LHP - YNCRight 0.015232640.06847509Inf 0.22245520.82395956
2LHP - YNCLeft -0.166647770.06847509Inf-2.43369920.01494541

Comparing acf between YNC & RHP¶

In [71]:
# SET REFERENCE
acfData2$Hand <- factor(acfData2$Hand, ordered = FALSE)
acfData2$Hand <- relevel(acfData2$Hand, "Right")
acfData2$Grp <- relevel(acfData2$Grp, "RHP")


# RUN ROBUST MODEL
robust.model2R = rlmer(acf_coeff ~ Hand + Grp + Hand:Grp + (1 | SubjID), acfData2 %>%
    filter(Grp != "LHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model2R, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model2R, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftRHP 0.25577350.07037451Inf 3.6344620.0002785613
2Right - LeftYNC-0.22643900.06094612Inf-3.7153970.0002028846
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1RHP - YNCRight 0.16113570.06582933Inf 2.4477791.437397e-02
2RHP - YNCLeft -0.32107690.06582933Inf-4.8774141.074859e-06

Relationship between trial-by-trial adaptive response and chronicity in RHP¶

In [72]:
fmchron2 = read.csv("fmchron_subj.csv", stringsAsFactors = FALSE)
acfData %>%
    filter(Grp != "YNC") %>%
    data.frame(., fmchron2[-22, ]) %>%
    filter(Grp == "RHP") %>%
    pivot_longer(cols = c(Left, Right), names_to = "Hand", values_to = "acf") %>%
    {
        forPlots2 <<- .
    }

colnames(forPlots2)[3] = "chron"
head(forPlots2)
A tibble: 6 × 7
SubjIDGrpchronuefmehiHandacf
<chr><chr><int><int><dbl><chr><dbl>
21RHP14159100Left -0.27089352
21RHP14159100Right 0.03707393
22RHP 9547 90Left -0.10022453
22RHP 9547 90Right 0.13830222
23RHP 6351100Left 0.04269213
23RHP 6351100Right 0.27618591
In [98]:
all = data.frame(acfData2,MMEstGains)
# all %>% filter(Grp!="LHP") %>% {.->>all2}
# head(all)
all$Hand = factor(all$Hand, levels=c('Left','Right'),labels = c("Left Hand","Right Hand"))

options(repr.plot.width=12, repr.plot.height=6.5)

Fig3a = ggplot(forPlots2, aes(y = acf, x = chron)) +
       geom_point(aes(col = Hand),stroke = 1.5, alpha = 0.9, size = 3.5, shape = 17, show.legend = FALSE) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       coord_capped_cart(xlim = c(0,200),ylim = c(-0.5,0.52),bottom = 'both', left = 'both') + 
       geom_smooth(aes(col=Hand,fill=Hand),  method = 'lm', formula = y~x, lwd=2,show.legend = FALSE,alpha=0.1) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + theme_pubclean() + thm +
       annotate_npc("paretic\n(right hand)", x = unit(c(0.85), "npc"), 
                    y = unit(c(0.61), "npc"), 
                    gp = gpar(fontsize = 18, fontface = 3,col = "#E93324")) +
       annotate_npc("non-paretic\n(left hand)", x = unit(c(0.85), "npc"), 
                    y = unit(c(0.18), "npc"),
                    gp = gpar(fontsize = 18, fontface = 3,col = "#110176")) +
       xlab("chronicity") + ylab("acf")

Fig3b = ggplot(all %>% filter(Grp!="LHP" & Hand=="Left Hand"), aes(x = acf_coeff, y = g)) +
       geom_point(aes(shape = Grp),stroke = 1.5, col="#110176",alpha = 1, size = 3.5, show.legend = FALSE) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_smooth(col="#110176",fill="#110176", method = 'lm', formula = y~x, lwd=2,show.legend = FALSE,alpha=0.1) +  
#        facet_grid(~Hand) +  #stat_ellipse(aes(group=Grp)) +
# stat_cor(method="spearman",size = 7) +
       coord_capped_cart(ylim = c(-0.8,0.2),xlim = c(-0.57,0.57), bottom = 'both', left = 'both') + 
      scale_shape_manual(values = c(16,17)) + 
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + 
   theme_pubclean() + thm +
       ylab("gain") + xlab("acf")

Fig4 = Fig3a + Fig3b + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 22))
Fig4
# ggsave(file="Fig4.svg", plot=Fig4, width=12, height=6.5)
In [75]:
chron_LRN = lmer(acf ~ chron:Hand + (1 | SubjID), data = forPlots2)
table4 = tab_model(chron_LRN, show.stat = TRUE)
display_html(head(table4$page.complete))
  acf
Predictors Estimates CI Statistic p
(Intercept) 0.20 0.03 – 0.37 2.33 0.020
chron * HandLeft -0.00 -0.00 – -0.00 -2.72 0.007
chron * HandRight 0.00 -0.00 – 0.00 0.19 0.849
Random Effects
σ2 0.02
τ00 SubjID 0.01
ICC 0.38
N SubjID 15
Observations 30
Marginal R2 / Conditional R2 0.385 / 0.617

Cursor error comparisons between groups¶

In [76]:
data  %>% 
filter(CondG==c("BmC")) %>% 
filter(Hand=="cursor")  %>%
select(Grp,Subj,CondG,SubNo,Cond,Dist,Hand,Trial,IDE,PVDE) %>% 
group_by(Subj,CondG) %>% 
mutate(meanIDE = mean(IDE),meanPVDE = mean(PVDE)) %>% 
distinct(Subj,Grp,CondG,meanIDE,meanPVDE) %>% 
{.->>cursorErrPlot}
cursorErrPlot$Grp = factor(cursorErrPlot$Grp, levels = c("LHP", "YNC", "RHP"))
cursorErrPlot$CondG <- as.factor(cursorErrPlot$CondG)
cursorErrPlot$CondG = factor(cursorErrPlot$CondG, levels = c("BmC", "BmS"),labels = c("One Cursor","Two Cursors"))

head(cursorErrPlot)
A grouped_df: 6 × 5
GrpSubjCondGmeanIDEmeanPVDE
<fct><chr><fct><dbl><dbl>
YNCSubj001One Cursor-3.1503930-0.3632224
YNCSubj002One Cursor-8.3391708-7.3970882
YNCSubj003One Cursor-1.1538569-0.6782084
YNCSubj004One Cursor-5.9896416-6.2597010
YNCSubj005One Cursor-0.2557512-1.1746078
YNCSubj006One Cursor-4.3837268-2.0762970
In [77]:
options(repr.plot.width=10, repr.plot.height=6.5)
comp = list(c("YNC","RHP"),c("YNC","LHP"),c("LHP","RHP"))

meanIDE = ggplot(data = cursorErrPlot,aes(x = Grp, y = meanIDE,shape=Grp)) +
geom_half_boxplot(size = 1,show.legend = FALSE) + geom_half_point_panel(size = 3,show.legend = FALSE) +
xlab("")+ ylab("cursor error (˚)") + 
coord_capped_cart(ylim = c(-20,20),bottom = "both", left = "both") +
stat_compare_means(comparisons = comp,label.y = c(12,15,18), na.rm = TRUE,label = "p.signif",
                   size=7,face="bold",bracket.size = 1,tip.length = 0.0) +  #facet_grid(~CondG)+
theme_pubclean() + thm

meanPVDE = ggplot(data = cursorErrPlot,aes(x = Grp, y = meanPVDE,shape=Grp)) +
geom_half_boxplot(show.legend = FALSE) + xlab("")+
coord_capped_cart(ylim = c(-12,12),bottom = "both", left = "both") +
stat_compare_means(comparisons = comp,label.y = c(5,7,9), na.rm = TRUE,label = "p.signif",
                   hide.ns = FALSE,vjust = 0, 
                    symnum.args=list(cutpoints = c(0.001, 0.01, 0.05,Inf),symbols = c("***", "**", "ns")),
                   size=8,face="bold",bracket.size = 1,tip.length = 0.0) +  
geom_half_point_panel(size = 3,show.legend = FALSE) + theme_pubclean() + thm
# ggsave(file="S3.svg", plot=meanIDE, width=10, height=6.5)


meanIDE | meanPVDE

Non-Paretic Left Hand Error¶

In [78]:
data  %>% 
filter(CondG==c("BmC")) %>% 
filter(Hand=="left")  %>%
filter(Grp!="LHP") %>% 
select(Grp,Subj,Hand,CondG,SubNo,Cond,Dist,Hand,Trial,IDE,PVDE) %>% 
group_by(Grp,Subj,Hand) %>%    
mutate(meanIDE = mean(IDE),meanPVDE = mean(PVDE)) %>% 
distinct(Subj,Grp,Hand,CondG,meanIDE,meanPVDE) %>% 
{.->>LhandErrPlot}
LhandErrPlot$Grp = factor(LhandErrPlot$Grp, levels = c("LHP", "YNC", "RHP"))

head(LhandErrPlot)
A grouped_df: 6 × 6
GrpSubjHandCondGmeanIDEmeanPVDE
<fct><chr><chr><chr><dbl><dbl>
YNCSubj001leftBmC -7.348676-11.3099125
YNCSubj002leftBmC-22.905913-21.5431340
YNCSubj003leftBmC 11.654869 2.6080192
YNCSubj004leftBmC -8.591172-15.7508540
YNCSubj005leftBmC 1.882509 -0.8127948
YNCSubj006leftBmC -5.524778 -9.9725552
In [79]:
options(repr.plot.width=8, repr.plot.height=6.5)
comp = list(c("YNC","RHP"))

LIDE = ggplot(data = LhandErrPlot,aes(x = Grp, y = meanIDE,shape=Grp)) +
geom_half_boxplot(col="#110176",fill="#110176",alpha = .3,size = 1,width = 0.6,show.legend = FALSE) + 
geom_half_point_panel(col="#110176",size = 3,show.legend = FALSE)+
geom_hline(yintercept = 0,col="black",lty="dashed",lwd=1) + 
xlab("")+ ylab("error at peak acc (˚)") + 
coord_capped_cart(ylim = c(-25,20),bottom = "both", left = "both") +
stat_compare_means(comparisons = comp,label.y = c(15), na.rm = TRUE,
                   size=7,face="bold",bracket.size = 1,tip.length = 0.0) + 
theme_pubclean() + thm

LPVDE = ggplot(data = LhandErrPlot,aes(x = Grp, y = meanPVDE,shape=Grp)) +
geom_half_boxplot(col="#110176",fill="#110176",alpha = 0.3,size = 1,width = 0.6,show.legend = FALSE) + 
geom_half_point_panel(col="#110176",size = 3,show.legend = FALSE)+
geom_hline(yintercept = 0,col="black",lty="dashed",lwd=1) + 
xlab("")+ ylab("error at peak speed (˚)") + 
coord_capped_cart(ylim = c(-25,20),bottom = "both", left = "both") +
stat_compare_means(comparisons = comp,label.y = c(15), na.rm = TRUE,
                   size=7,face="bold",bracket.size = 1,tip.length = 0.0) + 
theme_pubclean() + thm

LErrFig = LIDE | LPVDE

# ggsave(file="S4.svg", plot=LErrFig, width=8, height=6.5)

LErrFig

For PNAS Revisions¶


Within trial gains, direction of causality¶

In individual correlation plots, switching x and y (theta_cursor with delta_theta)¶

In [80]:
options(repr.plot.width = 13.5, repr.plot.height = 4.5)

Ind_subjGL = 
ggplot(bmCData2 %>% filter(Grp=="LHP" & Subj=="Subj006"), aes(y = Adjustment, x = IDE_cursor,na.rm = TRUE,group = Hand)) +
       geom_point(aes(col = Hand),stroke = 1.5, size = 3.5, shape = 15, show.legend = FALSE,alpha = 0.3) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'****'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),method="spearman",label.x = -20, label.y = c(20,15), size=8,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-20,20),bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand),  method = 'lm', fill=NA,formula = y~x, lwd=3.5,show.legend = FALSE,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + 
       scale_alpha_manual(values = c(0.7, 0.2)) +
       theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjL.svg', plot=Ind_subjGL, width=4.5, height=4.5)

Ind_subjG = 
ggplot(bmCData2 %>% filter(SubNo==14), aes(y = Adjustment, x = IDE_cursor,na.rm = TRUE)) +
       geom_point(aes(col = Hand),stroke = 1.5, alpha = 0.3, size = 3.5, shape = 16, show.legend = FALSE) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'***'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),label.x = -20, label.y = c(20,15),size=7.5,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-20,20), bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand,fill=Hand), fill=NA,method = 'lm', formula = y~x, lwd=3.5,show.legend = FALSE,alpha=0.1,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjG.svg', plot=Ind_subjG, width=4.5, height=4.5)

Ind_subjGR = 
ggplot(bmCData2 %>% filter(SubNo==22), aes(y = Adjustment, x = IDE_cursor,na.rm = TRUE,group = Hand)) +
       geom_point(aes(col = Hand),stroke = 1.5, size = 3.5, shape = 17, show.legend = FALSE,alpha = 0.3) + 
       geom_hline(yintercept = 0,size = 1,col="gray30",lty="dashed") + 
       geom_vline(xintercept = 0,size = 1,col="gray30",lty="dashed") + 
       stat_cor(aes(group = Hand,col=Hand,label =paste(..r.label.., cut(..p.., breaks = c(-Inf, 0.0001, 0.001, 0.01, 0.05, Inf),
                                                                        labels = c("'***'", "'***'", "'**'", "'*'", "'ns'")),
                                                       sep = "~")),digits=2,label.x = -12, label.y = c(20,15), size=7.5,show.legend = FALSE) + 
       coord_capped_cart(xlim = c(-20,20),ylim = c(-20,20),bottom = 'both', left = 'both') + #facet_wrap(~Hand,ncol=2) +
       geom_smooth(aes(col=Hand),  method = 'lm', fill=NA,formula = y~x, lwd=3.5,show.legend = FALSE,fullrange=TRUE) +  
       scale_color_manual(values=c("#110176","#E93324")) +
       scale_fill_manual(values=c("#110176","#E93324")) + 
       theme_pubclean() + thm +
       theme(panel.spacing.x=unit(-2, "lines")) + 
       theme(strip.text.x = element_text(size=0)) + theme(strip.background = element_blank(), strip.text = element_blank()) +
       ylab(expression(bold(Delta*theta))) + 
       xlab(expression(bold(theta)[Cursor]))

# ggsave(file='Ind_subjR.svg', plot=Ind_subjGR, width=4.5, height=4.5)

Ind_subjGALL = Ind_subjGL | Ind_subjG | Ind_subjGR

# ggsave(file='Ind_subjGALL.svg', plot=Ind_subjGALL, width=13.5, height=4.5)
Ind_subjGALL

In mixed effects models, switching x and y (adjustment with IDE_cursor).¶

  • Adjustment varies by Hand and moderates the relationship with IDE_cursor, so 'Hand' included as FE and RE.
  • Value of g does change from switching the x and y, but all main results still present.
  • Interestingly, it seems more clear that the negative gains are present for both hands in YNC, i.e., both hands contribute to correction within trial, just Left > Right. This is what was suggested by R2 as well & Dr. G.
  • In the right hemisphere stroke group, there is less compensation with both hands compared to YNC. This is why we refer to this as a "loss". Is there a better alternative word?
In [81]:
mixedModEstInter = suppressMessages(lmer(data = bmCData2, Adjustment ~   Hand + Grp + Grp:Hand:IDE_cursor +
    (1 + Hand + Hand:Grp:IDE_cursor| SubNo)))

coefs <- coef(summary(mixedModEstInter))
data.frame(coefs[, 1:2])
A data.frame: 10 × 2
EstimateStd..Error
<dbl><dbl>
(Intercept)-3.235790800.62795809
HandRight 5.616641720.86582681
GrpRHP-0.063778170.67888420
GrpLHP 0.060180700.70039911
HandLeft:GrpYNC:IDE_cursor-0.377246570.05020927
HandRight:GrpYNC:IDE_cursor-0.237201650.07155430
HandLeft:GrpRHP:IDE_cursor-0.099023400.04606562
HandRight:GrpRHP:IDE_cursor-0.407188470.07530568
HandLeft:GrpLHP:IDE_cursor-0.084741970.04442909
HandRight:GrpLHP:IDE_cursor-0.185751190.08265325
In [82]:
Left = c(coef(mixedModEstInter)$SubNo[1:20,"HandLeft:GrpYNC:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[21:32,"HandLeft:GrpRHP:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[44:47,"HandLeft:GrpRHP:IDE_cursor"],         
         coef(mixedModEstInter)$SubNo[33:43,"HandLeft:GrpLHP:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[48:50,"HandLeft:GrpLHP:IDE_cursor"])

Right = c(coef(mixedModEstInter)$SubNo[1:20,"HandRight:GrpYNC:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[21:32,"HandRight:GrpRHP:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[44:47,"HandRight:GrpRHP:IDE_cursor"],         
         coef(mixedModEstInter)$SubNo[33:43,"HandRight:GrpLHP:IDE_cursor"],
         coef(mixedModEstInter)$SubNo[48:50,"HandRight:GrpLHP:IDE_cursor"])


Grp = c(matrix('YNC',20),matrix('RHP',15),matrix('LHP',15)) #,

Subj = c(1:length(Left))
Slopes = data.frame(Subj,Grp,Left,Right)

# Pivot data table for analysis & plotting
Slopes  %>% 
pivot_longer(cols = c(Left,Right),names_to = "Hand",values_to = "g")  %>% 
{.->>MMEstGains}
MMEstGains$Grp = factor(MMEstGains$Grp, levels=c('LHP','YNC','RHP'))
head(MMEstGains)
A tibble: 6 × 4
SubjGrpHandg
<int><fct><chr><dbl>
1YNCLeft -0.45767043
1YNCRight-0.36023260
2YNCLeft -0.25721714
2YNCRight-0.06565796
3YNCLeft -0.74347454
3YNCRight-0.14888867
In [83]:
options(repr.plot.width = 11, repr.plot.height = 7)
options(warn = -1)

rs1 = ggplot(data = MMEstGains, aes(x = Grp, y = g, shape = Grp, col = Hand)) + stat_summary(mapping = NULL,
    data = NULL, fun.data = mean_se, alpha = 0.8, geom = "pointrange", position = position_dodge(width = 0.5),
    size = 1.5, lwd = 3, show.legend = FALSE) + geom_half_point_panel(aes(col = Hand),
    size = 3, alpha = 0.3, position = position_dodge(width = 0.5), show.legend = FALSE) +
    geom_line(aes(group = Hand, col = Hand), size = 2, fun.data = mean_se, alpha = 1,
        position = position_dodge(width = 0.5), stat = "summary", show.legend = FALSE) +
    geom_hline(yintercept = 0, linetype = "dashed", size = 1) + coord_capped_cart(ylim = c(-0.75,
    0.75), bottom = "both", left = "both") + xlab("") + ylab("gain") + annotate_npc("paretic left",
    x = unit(c(0.08), "npc"), y = unit(c(0.38), "npc"), rot = 90, gp = gpar(fontsize = 20,
        fontface = 3, col = "#110176")) + annotate_npc("paretic right", x = unit(c(0.92),
    "npc"), y = unit(c(0.28), "npc"), rot = 90, gp = gpar(fontsize = 20, fontface = 3,
    col = "#E93324")) + scale_color_manual(values = c("#110176", "#E93324")) + scale_fill_manual(values = c("#110176",
    "#E93324")) + scale_shape_manual(values = c(15, 16, 17)) + theme_pubclean() +
    thm + theme(axis.line.x = element_line(colour = "white"), axis.text.x = element_blank(),
    axis.ticks.x = element_blank(), legend.position = "none")


Fig1 = gr1 + rs1

# ggsave(file = "Fig1_gain.svg", plot = Fig1, width = 11, height = 7)

Fig1
options(warn = 0)

Robust LMER to compare g between groups and hands¶

Comparing g between YNC & LHP¶
In [84]:
# SET REFERENCE
MMEstGains$Hand <- factor(MMEstGains$Hand, ordered = FALSE)
MMEstGains$Hand <- relevel(MMEstGains$Hand, "Right")
MMEstGains$Grp <- relevel(MMEstGains$Grp, "LHP")

# RUN ROBUST MODEL
robust.model1L = rlmer(g ~ Hand + Grp + Hand:Grp + (1 | Subj), MMEstGains %>% filter(Grp != "RHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model1L, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model1L, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftLHP-0.13514330.04148137Inf-3.2579291.122285e-03
2Right - LeftYNC 0.15283620.03592392Inf 4.2544412.095717e-05
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1LHP - YNCRight0.0066369810.05703829Inf0.11636019.073672e-01
2LHP - YNCLeft 0.2946165120.05703829Inf5.16524102.401294e-07
Comparing g between YNC & RHP¶
In [85]:
# SET REFERENCE
MMEstGains$Hand <- factor(MMEstGains$Hand, ordered = FALSE)
MMEstGains$Hand <- relevel(MMEstGains$Hand, "Right")
MMEstGains$Grp <- relevel(MMEstGains$Grp, "RHP")

# RUN ROBUST MODEL
robust.model1R = rlmer(g ~ Hand + Grp + Hand:Grp + (1 | Subj), MMEstGains %>% filter(Grp != "LHP"), doFit = TRUE)

#EMMEANS
summary(emmeans(robust.model1R, specs = pairwise ~ Hand | Grp)$contrasts)
summary(emmeans(robust.model1R, specs = pairwise ~ Grp | Hand)$contrasts)
A summary_emm: 2 × 7
contrastGrpestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1Right - LeftRHP-0.28465530.04341018Inf-6.5573395.477632e-11
2Right - LeftYNC 0.15455370.03759432Inf 4.1110923.937915e-05
A summary_emm: 2 × 7
contrastHandestimateSEdfz.ratiop.value
<fct><fct><dbl><dbl><dbl><dbl><dbl>
1RHP - YNCRight-0.13669320.05097017Inf-2.6818287.322102e-03
2RHP - YNCLeft 0.30251580.05097017Inf 5.9351532.935711e-09

Left hand errors result from an inability to compensate for inertial properties and interjoint dynamics as shown in previous work. Our data in neurotypical controls support this.¶

In [86]:
data %>% 
filter(CondG==c("BmC","BmS")) %>% 
filter(Hand!="cursor") %>% 
filter(Grp=="YNC") %>% 
mutate(Trial = as.numeric(Trial)) %>% 
select(Subj,SubNo,CondG,Dist,Hand,Trial,IDE,PVDE,XPosPV,YPosPV,CURV) %>% 
group_by(Subj,CondG,Hand) %>% 
summarize(mnIDE = median(IDE),
          mnPVDE = median(PVDE),
          mnCurv = median(CURV),.groups="drop") %>% 
{.->>LErrorsDat}

head(LErrorsDat)

LErrorsDat$CondG = factor(LErrorsDat$CondG, levels = c("BmC", "BmS"),labels = c("One Cursor","Two Cursors"))
A tibble: 6 × 6
SubjCondGHandmnIDEmnPVDEmnCurv
<chr><chr><chr><dbl><dbl><dbl>
Subj001BmCleft -6.678415-10.4204800.169935
Subj001BmCright -0.110695 4.7595700.126945
Subj001BmSleft -6.329610-10.8530100.169910
Subj001BmSright -4.097490 3.8341250.119385
Subj002BmCleft -23.366285-22.3714950.294570
Subj002BmCright 10.431790 10.0185700.123570

Initial directional errors and curvature of the left hand from our neurotypical data¶

In [87]:
comp = list(c("left","right"))

theta_peakAcc = 
ggplot(LErrorsDat, aes(y=(mnIDE),x=Hand,col=Hand)) + 
  geom_half_dotplot(aes(col=Hand,fill=Hand),binwidth=1.5,alpha = 0.7,show.legend = FALSE)+ 
  geom_half_boxplot(aes(col=Hand,fill=Hand),alpha = 0.2,size = 1, show.legend = FALSE)+ 
coord_capped_cart(ylim= c(-30, 20),bottom = 'both', left='both') + facet_grid(~CondG) +
geom_hline(yintercept = 0,lwd=1.5,lty="solid",col="black") + ylab("initial directional error") + xlab("") +
scale_color_manual(values = c("#110176", "#E93324")) + scale_fill_manual(values = c("#110176","#E93324")) +
stat_compare_means(comparisons = comp,label.y = c(15), na.rm = TRUE,
                   size=7,face="bold",bracket.size = 1,tip.length = 0.0) + 
theme_pubclean() + theme(axis.line.x = element_line(colour = "white"), axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
 ylab(expression(bold(theta)[italic(peak~acc)])) + ggtitle("A. Directional error is larger for the left hand") +
thm

curvature = 
ggplot(LErrorsDat, aes(y=(mnCurv),x=Hand,col=Hand)) + 
  geom_half_dotplot(aes(col=Hand,fill=Hand),binwidth=0.015,alpha = 0.7,show.legend = FALSE)+ 
  geom_half_boxplot(aes(col=Hand,fill=Hand),alpha = 0.2,size = 1, show.legend = FALSE)+ 
coord_capped_cart(ylim= c(0, 0.4),bottom = 'both', left='both') + facet_grid(~CondG) +
scale_color_manual(values = c("#110176", "#E93324")) + scale_fill_manual(values = c("#110176","#E93324")) +
stat_compare_means(comparisons = comp,label.y = c(0.34), na.rm = TRUE,
                   size=7,face="bold",bracket.size = 1,tip.length = 0.0) + 
theme_pubclean() + theme(axis.line.x = element_line(colour = "white"), axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
 ylab("path curvature (au)") + xlab("") + ggtitle("B. Hand path is more curved for the left hand") + thm

options(repr.plot.width=8, repr.plot.height=11)
newSupp1 = theta_peakAcc / curvature

# ggsave(file = "S2.svg", plot = newSupp1, width = 8, height = 11)
newSupp1

Time between moments of peak acc and peak speed and its comparison between groups¶

In [88]:
# colnames(data)[1] ="Grp"
data$Grp = factor(data$Grp, levels=c('YNC','RHP','LHP'))


data %>% 
filter(CondG==c("BmC")) %>% 
mutate(Adj = PVDE - IDE,FFDur = (TPV-PPT)*1000) %>%
filter(Hand!="cursor") %>% 
mutate(Trial = as.numeric(Trial)) %>% 
select(Grp,SubNo,Cond,Dist,Hand,Trial,IDE,PPT,Adj,FFDur) %>% 
group_by(Grp,SubNo,Hand) %>% 
summarize(mnAdj = median(Adj),
          mnFFDur = median(FFDur),.groups="drop") %>% 
{.->>AdjFF}

head(AdjFF)

AdjFF %>% 
group_by(Grp) %>% 
summarize(GrpMeanFF = median(mnFFDur), GrpSEFF = parameters::standard_error(mnFFDur))
A tibble: 6 × 5
GrpSubNoHandmnAdjmnFFDur
<fct><int><chr><dbl><dbl>
YNC1left -4.33657051.72138
YNC1right 6.06369560.34000
YNC2left 1.28054551.72069
YNC2right-0.33088060.34000
YNC3left -6.81613043.10138
YNC3right 7.16192043.10069
A tibble: 3 × 3
GrpGrpMeanFFGrpSEFF
<fct><dbl><dbl>
YNC 51.72069 1.547207
RHP131.4662124.583954
LHP142.2427622.948625
In [89]:
# Define function to calculate IQR at given quantiles
iqr = function(z, lower = 0.25, upper = 0.75) {
  data.frame(
    y = median(z),
    ymin = quantile(z, lower),
    ymax = quantile(z, upper)
  )
}
In [90]:
iqr(AdjFF$mnFFDur[AdjFF$Grp=="RHP"])
range(AdjFF$mnFFDur[AdjFF$Grp=="RHP"])
A data.frame: 1 × 3
yyminymax
<dbl><dbl><dbl>
25%131.466261.42573167.0239
  1. -301.724138
  2. 254.3127585
In [91]:
options(repr.plot.width=6, repr.plot.height=6)
my_comparisons = list(c("YNC","LHP"),c("YNC","RHP"),c("LHP","RHP"))
options(warn=-1)
AdjFFDurFig = 
ggplot(aes(x = Grp,y =(mnFFDur),shape=Grp),data =AdjFF) + 
stat_summary(fill="gray",geom="bar",width=.75, alpha = 0.8,fun = median,show.legend = FALSE) + 
stat_summary(fun.data = iqr, fun.args = list(lower = 0.25, upper = 0.75),  lwd=2, geom = "errorbar",width=0.2, show.legend = FALSE) + 
geom_point(aes(col=Hand,shape=Grp),na.rm=TRUE, alpha = 0.5,position = position_jitterdodge(jitter.width = 0.3,jitter.height = 0.2,dodge.width = 0.75),size=4,show.legend=FALSE) +
coord_capped_cart(ylim= c(0, 290),bottom = 'both', left='both') + 
scale_color_manual(values = c("#110176", "#E93324")) + scale_fill_manual(values = c("#110176","#E93324")) +
xlab("") + ylab("Early Phase\nDuration (ms)") +
ylab(expression(t[italic(peak~speed)]-t[italic(peak~acc)]~(ms))) +
stat_compare_means(show.legend=FALSE,vjust = 0.5,comparisons = my_comparisons,label.y=c(200,245,270), na.rm = TRUE,
                   size=12,face="italic",bracket.size = 1,tip.length = 0.0, hide.ns = TRUE,
                  symnum.args=list(cutpoints = c(0.0001, 0.001, 0.05,Inf),symbols = c("***", "**", "ns"))) +
theme_pubclean() + thm + theme(legend.position="none")
AdjFFDurFig
# ggsave(file = "AdjFFDurFig.svg", plot = AdjFFDurFig, width = 6, height = 6)

Figure 1E & F (+Supp Stroke): one-cursor vs two-cursor mean-centered in YNC, LHP and RHP¶

In [92]:
incl = colnames(data)[- (1:8)]
data  %>% 
filter(CondG==c("BmC","BmS")) %>% 

unite(col = "group", Grp,Subj,Cond,Trial) %>% 
filter(!group %in% rem2) %>% 
separate(col = group, into = c("Grp","Subj","Cond","Trial"),sep = "_") %>% 
pivot_wider(id_cols = NULL, names_from= Hand,values_from = all_of(incl)) %>% 
{.->>newData}
tail(newData,5)

newData$Grp = factor(newData$Grp, levels=c('YNC','RHP','LHP')) 
colnames(newData)
A tibble: 5 × 79
GrpSubjCondTrialSubNoCondGDistCT_leftCT_rightCT_cursor⋯YPosPV_cursorXPosOff_leftXPosOff_rightXPosOff_cursorYPosOff_leftYPosOff_rightYPosOff_cursorXDistPV_leftXDistPV_rightXDistPV_cursor
<chr><chr><chr><chr><int><chr><int><dbl><dbl><dbl>⋯<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
LHPSubj016BmC151151BmC150.9230.9230.923⋯1.1194680.92562641.0716260.99932871.0745701.1569621.127890-0.019821821-0.013992917-0.08600599
LHPSubj016BmC151351BmC150.8850.8850.885⋯1.1051350.93468891.0720091.00275851.1242451.1196741.114799-0.013083113-0.015627742-0.08395930
LHPSubj016BmC151551BmC150.9950.9950.995⋯1.0807340.92733291.0659190.99655271.0853241.0990061.092654-0.016961871-0.020339184-0.08830897
LHPSubj016BmC151751BmC150.9280.9280.928⋯1.0676730.94623091.0814731.01179211.0782811.0973191.074236-0.001889574-0.005672671-0.07449025
LHPSubj016BmC151951BmC150.9230.9230.923⋯1.1112230.93388251.0688111.00097891.1215911.1249811.128809-0.008924704-0.016251533-0.08224652
  1. 'Grp'
  2. 'Subj'
  3. 'Cond'
  4. 'Trial'
  5. 'SubNo'
  6. 'CondG'
  7. 'Dist'
  8. 'CT_left'
  9. 'CT_right'
  10. 'CT_cursor'
  11. 'RT_left'
  12. 'RT_right'
  13. 'RT_cursor'
  14. 'PPT_left'
  15. 'PPT_right'
  16. 'PPT_cursor'
  17. 'TPV_left'
  18. 'TPV_right'
  19. 'TPV_cursor'
  20. 'MT_left'
  21. 'MT_right'
  22. 'MT_cursor'
  23. 'PV_left'
  24. 'PV_right'
  25. 'PV_cursor'
  26. 'EE_left'
  27. 'EE_right'
  28. 'EE_cursor'
  29. 'VE_left'
  30. 'VE_right'
  31. 'VE_cursor'
  32. 'IDE_left'
  33. 'IDE_right'
  34. 'IDE_cursor'
  35. 'PRDE_left'
  36. 'PRDE_right'
  37. 'PRDE_cursor'
  38. 'PVDE_left'
  39. 'PVDE_right'
  40. 'PVDE_cursor'
  41. 'FDE_left'
  42. 'FDE_right'
  43. 'FDE_cursor'
  44. 'FPE_left'
  45. 'FPE_right'
  46. 'FPE_cursor'
  47. 'IDIST_left'
  48. 'IDIST_right'
  49. 'IDIST_cursor'
  50. 'PVDIST_left'
  51. 'PVDIST_right'
  52. 'PVDIST_cursor'
  53. 'FDIST_left'
  54. 'FDIST_right'
  55. 'FDIST_cursor'
  56. 'CURV_left'
  57. 'CURV_right'
  58. 'CURV_cursor'
  59. 'OnAxErr_left'
  60. 'OnAxErr_right'
  61. 'OnAxErr_cursor'
  62. 'OffAxErr_left'
  63. 'OffAxErr_right'
  64. 'OffAxErr_cursor'
  65. 'XPosPV_left'
  66. 'XPosPV_right'
  67. 'XPosPV_cursor'
  68. 'YPosPV_left'
  69. 'YPosPV_right'
  70. 'YPosPV_cursor'
  71. 'XPosOff_left'
  72. 'XPosOff_right'
  73. 'XPosOff_cursor'
  74. 'YPosOff_left'
  75. 'YPosOff_right'
  76. 'YPosOff_cursor'
  77. 'XDistPV_left'
  78. 'XDistPV_right'
  79. 'XDistPV_cursor'
In [93]:
## FOR SPATIAL COOUPLING FIGURE

rXPos_bmC = vector()
rXPos_bmS = vector()

SubjID = vector()

for (i in newData$SubNo)
    {rXPos_bmC[i] = cor.test(newData$XDistPV_left[newData$SubNo==i & newData$CondG=="BmC"],
                             newData$XDistPV_right[newData$SubNo==i & newData$CondG=="BmC"],
                             method = "spearman")$ estimate
     rXPos_bmS[i] = cor.test(newData$XDistPV_left[newData$SubNo==i & newData$CondG=="BmS"],
                             newData$XDistPV_right[newData$SubNo==i & newData$CondG=="BmS"],
                             method = "spearman")$ estimate

     SubjID[i] = toString(i)
    }


rXPos = data.frame(rXPos_bmC,rXPos_bmS)


colnames(rXPos)[1] <- "BmC"
colnames(rXPos)[2] <- "BmS"


j = c(matrix('YNC',20),matrix('RHP',15),matrix('LHP', 16)) #,

rData = data.frame(SubjID,data.frame(j), rXPos)

colnames(rData)[2] = "Grp"

rData  %>% 
pivot_longer(cols = c("BmC","BmS"),"Cond",values_to = "corr")  %>% 

{.->> rData2}

rData2$Grp = factor(rData2$Grp, levels=c('YNC','RHP','LHP'))

head(rData2,4)
A tibble: 4 × 4
SubjIDGrpCondcorr
<chr><fct><chr><dbl>
1YNCBmC-0.9157895
1YNCBmS-0.6571429
2YNCBmC-0.9639098
2YNCBmS-0.3669173
In [94]:
display_markdown('### EXP vs CTRL condition in YNC')
options(repr.plot.width=12, repr.plot.height=6)

glance(t.test(rData2$corr[rData2$Grp=="YNC"&rData2$Cond=="BmS"],
       rData2$corr[rData2$Grp=="YNC"&rData2$Cond=="BmC"],paired = TRUE))
C1 = 
ggplot(aes(x = (XDistPV_right)*100, y=(XDistPV_left)*100, col=CondG),
           data = newData %>% filter(Grp=="YNC")) + 
    geom_point(size=2.5,alpha = 0.7) + stat_ellipse(lwd=2) + 
    coord_capped_cart(xlim = c(-6,6),ylim = c(-6,6), bottom = 'both', left='both') + 
    geom_abline(slope = -1,lwd=1) + 
    ylab("") + xlab("") + 
    theme_pubclean() + thm +  facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"))

D1 = 
ggplot(rData2 %>% filter(Grp=="YNC"), aes(x=corr)) + 
    geom_density(aes(y=..scaled..,colour=Cond, fill=Cond),alpha = 0.5, size=2,show.legend = TRUE) + 
    geom_rug(aes(col = Cond),size = 1.4,length = unit(0.07, "npc"),outside = FALSE, sides = "b") + 
    geom_vline(xintercept = 0,col="black",size=1,linetype = "dashed")+ 
    coord_capped_cart(ylim = c(-0.1,1), xlim = c(-1, 1), bottom = 'both', left='both') + 
    ylab("density") + theme_pubclean() + thm + facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"),name="", labels = c("One cursor","Two cursors"))

YNC = C1 + D1 + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 22))
YNC
# ggsave(file = "YNC_1_2_corr.svg", plot = YNC, width = 12, height = 12)

EXP vs CTRL condition in YNC¶

A tibble: 1 × 8
estimatestatisticp.valueparameterconf.lowconf.highmethodalternative
<dbl><dbl><dbl><dbl><dbl><dbl><chr><chr>
0.46406026.3705114.125443e-06190.31159370.6165266Paired t-testtwo.sided
In [95]:
display_markdown('## EXP vs CTRL condition in RHP')
options(repr.plot.width=12, repr.plot.height=6)

glance(t.test(rData2$corr[rData2$Grp=="RHP"&rData2$Cond=="BmS"],
       rData2$corr[rData2$Grp=="RHP"&rData2$Cond=="BmC"],paired = TRUE))

C2 = 
ggplot(aes(x = (XDistPV_right)*100, y=(XDistPV_left)*100, col=CondG),
           data = newData %>% filter(Grp=="RHP")) + 
    geom_point(size=2.5,alpha = 0.7) + stat_ellipse(lwd=2) + 
    coord_capped_cart(xlim = c(-6,6),ylim = c(-6,6), bottom = 'both', left='both') + 
    geom_abline(slope = -1,lwd=1) + 
    ylab("") + xlab("") + 
    theme_pubclean() + thm +  facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"))

D2 = 
ggplot(rData2 %>% filter(Grp=="RHP"), aes(x=corr)) + 
    geom_density(aes(y=..scaled..,colour=Cond, fill=Cond),alpha = 0.5, size=2,show.legend = TRUE) + 
    geom_rug(aes(col = Cond),size = 1.4,length = unit(0.07, "npc"),outside = FALSE, sides = "b") + 
    geom_vline(xintercept = 0,col="black",size=1,linetype = "dashed")+ 
    coord_capped_cart(ylim = c(-0.1,1), xlim = c(-1, 1), bottom = 'both', left='both') + 
    ylab("density") + theme_pubclean() + thm + facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"),name="", labels = c("One cursor","Two cursors"))

RHP = C2 + D2 + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 22))
RHP

EXP vs CTRL condition in RHP¶

A tibble: 1 × 8
estimatestatisticp.valueparameterconf.lowconf.highmethodalternative
<dbl><dbl><dbl><dbl><dbl><dbl><chr><chr>
0.55044284.4211120.0005807458140.28340990.8174757Paired t-testtwo.sided
In [96]:
display_markdown('## EXP vs CTRL condition in LHP')
options(repr.plot.width=12, repr.plot.height=6)

glance(t.test(rData2$corr[rData2$Grp=="LHP"&rData2$Cond=="BmS"],
       rData2$corr[rData2$Grp=="LHP"&rData2$Cond=="BmC"],paired = TRUE))

C3 = 
ggplot(aes(x = (XDistPV_right)*100, y=(XDistPV_left)*100, col=CondG),
           data = newData %>% filter(Grp=="LHP")) + 
    geom_point(size=2.5,alpha = 0.7) + stat_ellipse(lwd=2) + 
    coord_capped_cart(xlim = c(-6,6),ylim = c(-6,6), bottom = 'both', left='both') + 
    geom_abline(slope = -1,lwd=1) + 
    ylab("") + xlab("") + 
    theme_pubclean() + thm +  facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"))

D3 = 
ggplot(rData2 %>% filter(Grp=="LHP"), aes(x=corr)) + 
    geom_density(aes(y=..scaled..,colour=Cond, fill=Cond),alpha = 0.5, size=2,show.legend = TRUE) + 
    geom_rug(aes(col = Cond),size = 1.4,length = unit(0.07, "npc"),outside = FALSE, sides = "b") + 
    geom_vline(xintercept = 0,col="black",size=1,linetype = "dashed")+ 
    coord_capped_cart(ylim = c(-0.1,1), xlim = c(-1, 1), bottom = 'both', left='both') + 
    ylab("density") + theme_pubclean() + thm + facet_wrap(~Grp) +
    scale_color_manual(values=c("#503010","#B5833F"),name = "",labels = c("One cursor","Two cursors")) + 
    scale_fill_manual(values=c("#503010","#B5833F"),name="", labels = c("One cursor","Two cursors"))

LHP = C3 + D3 + plot_annotation(tag_levels = 'A') & theme(plot.tag = element_text(size = 22))
LHP

EXP vs CTRL condition in LHP¶

A tibble: 1 × 8
estimatestatisticp.valueparameterconf.lowconf.highmethodalternative
<dbl><dbl><dbl><dbl><dbl><dbl><chr><chr>
0.58372615.4801368.106093e-05140.35527040.8121817Paired t-testtwo.sided

END OF ANALYSIS¶


In [97]:
display_markdown("### Packages in use: ")
search()

display_markdown("### R version & citation: ")
citation()

display_markdown("### Specific package to cite: EMMEANS")
citation("emmeans")

Packages in use:¶

  1. '.GlobalEnv'
  2. 'package:ggplotify'
  3. 'package:influence.ME'
  4. 'package:svglite'
  5. 'package:broom'
  6. 'package:emmeans'
  7. 'package:table1'
  8. 'package:sjstats'
  9. 'package:sjPlot'
  10. 'package:nlme'
  11. 'package:grid'
  12. 'package:ggExtra'
  13. 'package:robustlmm'
  14. 'package:gghalves'
  15. 'package:lme4'
  16. 'package:Matrix'
  17. 'package:lemon'
  18. 'package:tidyr'
  19. 'package:dplyr'
  20. 'package:patchwork'
  21. 'package:ggpubr'
  22. 'package:ggplot2'
  23. 'package:repr'
  24. 'package:IRdisplay'
  25. 'package:jsonlite'
  26. 'package:formatR'
  27. 'jupyter:irkernel'
  28. 'package:stats'
  29. 'package:graphics'
  30. 'package:grDevices'
  31. 'package:utils'
  32. 'package:datasets'
  33. 'package:methods'
  34. 'Autoloads'
  35. 'package:base'

R version & citation:¶

To cite R in publications use:

  R Core Team (2020). R: A language and environment for statistical
  computing. R Foundation for Statistical Computing, Vienna, Austria.
  URL https://www.R-project.org/.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {R: A Language and Environment for Statistical Computing},
    author = {{R Core Team}},
    organization = {R Foundation for Statistical Computing},
    address = {Vienna, Austria},
    year = {2020},
    url = {https://www.R-project.org/},
  }

We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also ‘citation("pkgname")’ for
citing R packages.

Specific package to cite: EMMEANS¶

To cite package ‘emmeans’ in publications use:

  Russell V. Lenth (2021). emmeans: Estimated Marginal Means, aka
  Least-Squares Means. R package version 1.7.0.
  https://CRAN.R-project.org/package=emmeans

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
    author = {Russell V. Lenth},
    year = {2021},
    note = {R package version 1.7.0},
    url = {https://CRAN.R-project.org/package=emmeans},
  }
In [ ]: